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SECTION  I 


INTRODUCTION 


Self  excited  oscillations  have  been  experimentally  observed  in  sepa¬ 
rated  flows  for  over  hundred  years.  Rayleigh  [1]  in  1880  proved  that  for 
inviscid,  incompressible  flow  the  unstable  velocity  profiles  must  have  an 
inflection  point.  Tollmein  [2]  in  1935  showed  that  for  symmetrical  velo¬ 
city  distributions,  or  for  velocity  distributions  of  the  boundary  layer 
type,  the  existence  of  the  inflection  point  implies  instability. 

Recently  Hankey  and  Shang  [3]  have  examined  the  self  induced  pressure 
oscillations  in  an  open  cavity.  Their  numerical  computations  compare  very 
well  with  the  previous  experimental  investigations.  Roscoe  and  Hankey  [4] 
have  studied  the  stability  of  hyperbolic  tangent  velocity  profile  in  a 
compressible  fluid,  while  Hankey,  Hunter  and  Harney  [5]  have  examined  the 
self-sustained  oscillations  (Buzz)  on  spiked  tipped  bodies  for  large  Mach 
numbers.  However,  a  systematic  stability  analysis  of  separated  flows  has 
not  been  undertaken.  It  is  the  purpose  of  this  report  to  conduct  a  sta¬ 
bility  analysis  of  a  general  class  of  separated  flows  (i.e.  reversed  flow 
Falkner-Skan)  in  order  to  help  shed  light  on  the  phenomenon  of  self-excited 
fluid  flows. 
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SECTION  II 


OBJECTIVES  OF  THE  RESEARCH  EFFORT 

The  objective  of  this  research  effort  was  to  analyze  a  series  of 
similar  separated  flows  for  different  values  of  3  and  to  determine  the 
amplification  factors  and  propagation  velocities  in  all  these  different 
cases.  Eight  cases  of  different  3  were  identified  to  be  analyzed.  These 
cases  were  those  with  reversed  flows  which  contained  velocity  profiles 


with  inflection  points. 


SECTION  III 


MEAN  :7L0M  EQUATIONS 


In  this  report,  incompressible  flows  will  be  analyzed.  In  subsequent 
work  we  plan  to  analyze  the  compressible  flows. 

The  incompressible  two-dimensional  Navier-Stokes  equations  are  as 
follows 

1  o 

(3.1) 

(3.2) 

(3.3) 

x  y 

Applying  the  boundary  layer  approximations  to  the  above  equations  for 
steady  flows  results  in  the  following: 

(3. A) 


u  + 

UU  +  VU  = 

_  1 

p 

+  vV2U 

t 

x  y 

P 

X 

V  + 

UV  +  W  = 

1 

P 

+  vV2V 

t 

x  y 

”  P 

y 

U  +  V  = 

0 

UU  +  VU  =  U  U  +  VU 
x  y  e  ex  yy 


U  +  V  =  0 
x  y 


(3.5) 


These  equations  may  be  reduced  to  one  ordinary  differential  equation  for 


the  case  where  Ug  =  cx  by  transforming  with  similarity  variables. 


d£  = 


Uedx 


(3.6) 


Hence 


where 


and 


dn  =  - 


Uedy 


/2^  v 


:,2 


f"'  +  ff"  -  0 (f  * “  -  1) 

u 


f'(n)  = 


u 


A  2m 

B  "  0  ^TTT  =  constant 

e 


with  boundary  conditions 

f(0)  -  0,  f '  (0)  -  0,  f  •  (oo)  -  i 


(3.7) 

(3.8) 

(3.9) 

(3.10) 

(3.11) 
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Falkner  and  Skan  [6]  originally  derived  this  equation  for  attached  flows 
however,  Stewartson  [7]  discovered  a  lower  branch  to  these  solutions  which 
represented  reversed  flows  from  incipient  separation  to  the  Chapman  solu¬ 
tion.  Christian,  Hankey  and  Petty  [8]  have  tabulated  these  solutions  for 
compressible  and  incompressible  flows.  It  is  this  wide  class  of  flows 
(which  have  inflection  points)  that  are  known  to  be  unstable  for  which 
we  shall  now  perform  a  stability  analysis. 
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SECTION  IV 


PERTURBATION  EQUATIONS 


Let  us  assume  small  perturbations  of  the  form 

ia(x  -  ct) 


U  =  U  (y)  +  U(y)  e 


V  = 


$(y)  e 


ia(x  -  ct) 


p  =  Pe(x)  +  P(y)  e 


ia(x  -  ct) 


(4.1) 

(4.2) 

(4.3) 


where  c  =  c  +  ic^  and  U,  <()  and  P  are  small  in  comparison  to  the  mean 
quantities.  If  we  substitute  these  values  of  U,  V,  and  P  in  equations 
(3.1),  (3.2)  and  (3.3);  retain  only  the  first  order  terms  and  assume  that 
Reynolds  number  Uex/v  is  large  then  the  equations  (4.1),  (4.2)  and  (4.3) 


reduce  to  one  single  equation 


>"  -  (ct  + 


U" 


■)♦  =  0 


(4.4) 


U  -  c 

The  classical  Rayleigh  equation  with  the  boundary  conditions 

<f>(0)  =  0,  <f>(°°)  =  0  (4 . 5a, b) 

By  transforming  the  equation  from  y  to  the  0  variable  we  obtain  the  fol¬ 
lowing  equation 

<f>nn  '  (“2  +  f~f-'c  )<f>  =  0  (4,6) 

where 


a  =  a  ^ 
dr) 

By  inserting  the  values  of  f ' (ri ,  3)  into  the  Rayleigh  equation  c(a,  3) 
can  be  determined  as  an  eigenvalue  which  satisfies  the  boundary  conditions 
(4.5a,b) . 
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SECTION  V 


SOLVING  SCHEME 


Eigenvalues  were  determined  by  a  shooting  method;  starting  with  a 
given  boundary  condition,  integrating  over  the  range  of  n  and  comparing 
the  result  with  the  outer  boundary  condition,  namely  i  =  0  at  n 

max 

The  process  involved  minimization  of  the  error  in  the  outer  boundary 

condition  which  was  chosen  to  be  the  square  of  the  norm  of  $ , 

2  2  2 

:•  =  0  +  .  =  SSQ.  (See  Appendix  3).  The  integration  was  done  using 

a  fourth-order  Runge-Kutta  method. 

The  method  of  finding  eigenvalues  utilized  a  minimization  routine 
written  primarily  by  Roscoe  [4],  Starting  from  a  given  guess  the 
routine  searched  along  a  constant  line  of  c^  with  increasing  steps  until 
it  found  a  relative  minimum  of  the  error.  It  then  used  the  last  three 
calculated  points  to  determine  a  parabola,  with  the  c  value  at  the  vertex 
used  as  the  latest  approximation.  Then  this  value  of  cr  was  held  constant 
and  a  search  along  a  line  of  changing  c^  was  carried  out.  After  a  new 
minimum  was  found,  the  quadratic  approximation  was  again  used  to  determine 
a  new  value  for  c^.  The  third  step  involved  searching  the  line  connecting 
the  original  guess  and  the  new  point.  After  finding  a  minimum  and  utilizing 
quadratic  approximation,  the  error  was  checked  to  see  if  it  was  less  than 
some  preset  limit.  If  not,  the  routine  started  again  with  the  latest  value 
used  in  place  of  the  original  guess. 

Generally,  the  routine  worked  quite  well.  Most  of  the  search  time 
was  attributable  to  bad  guesses  and  finding  the  direction  in  which  the 
search  should  be  continued.  An  eigenvalue  was  usually  located  in  a  very 
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narrow  region  of  the  plane  and  even  though  the  step  size  was  continually 

reduced,  it  was  frequently  large  enough  to  move  the  test  point  out  of  the 

acceptable  region.  For  example,  the  initial  guess  in  one  case  led  to  an 
13 

error  of  4.1  x  10  ,  however,  after  only  128  new  error  computations,  the 

-4 

error  had  been  reduced  17  orders  of  magnitude  to  1.9  x  10  ,  while  cr  had 

been  changed  by  4.25%  and  c^  had  been  changed  by  3.82%.  Convergence  was 
also  retarded  for  small  values  of  c^,  e.g.  jc_J  <  .001.  This  was  concur¬ 
rent  with  c^  approaching  its  limiting  value. 

The  Howard  semicircle  theorem  [9]  was  used  as  an  aid  in  determining 
suitable  initial  guesses.  If  c^  is  the  propagation  velocity,  a  is  the 

wave  number,  c.  is  the  amplification  factor,  and  U  and  U  .  are  the 
i  max  min 

maximum  and  minimum  values  of  the  range  of  U,  the  theorem  states 

[c  -  1/2(U  +  U  .  )]2  +  c2.  <  [1/2(U  -  U  .  )12. 

r  max  min  l  —  '  max  min' 1 

Thus,  the  complex  wave  velocity  for  an  unstable  mode  lies  inside  the  upper 

semi-circle  which  has  the  range  of  U  as  diameter. 


SECTION  VI 


RESULTS 


Eight  cases  were  computed  for  8  values  of  -.0001,  -.0005,  -.002, 
-.04,  -.08,  -.12,  -.16  and  -.19884.  For  a  wide  range  of  a  values  the 
eigenvalues  were  ascertained.  These  values  are  tabulated  in  tables  B-l  - 
B-8  in  Appendix  B.  a  is  related  to  a  by  the  relation 


* 

a6 


a  (1  -  f')  dn 

o 


* 


The  values  of  and  versus  a6  are  plotted 
Figure  3a-3h  shows  Howard's  plot  (9)  for  these 
eigenvalues  for  a  series  of  solutions  are  also 
Appendix  B. 


in  figures  la-lh  and  2a-2h. 
solutions.  Some  typical 
tabulated  and  plotted  in 
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LEGEND 


LEGEND 

LETTER  /3 

a 

-.0001 

b 

-.0005 

c 

-.002 

d 

-.04 

e 

-.08 

f 

-.12 

i 

-.15 

h 

-.19884 

SECTION  VII 


CONCLUSIONS 


The  stability  of  a  series  of  similar  separated  flows  have  been 
analyzed.  Amplification  factors  and  propagation  velocities  of  the  dis¬ 
turbances  were  determined.  The  results  show  that  a  small  zone  of  insta¬ 
bility  does  exist  for  these  flows  with  inflexion  points.  The  ampli- 
faction  factor  increases  as  the  extent  of  the  reversed  flow  increases. 
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SECTION  VIII 


RECOMMENDATIONS 


Suggestions  for  follow-on  research:  We  would  like  to  investigate 
the  instability  of  laminar  separated  flows  under  the  influence  of  compres¬ 
sibility.  For  a  hyperbolic  tangent  velocity  profile  Roscoe  (4)  showed 
the  instability  to  diminish  with  the  increase  of  Mach  number  until  the 
Rayleigh  instability  actually  vanished  at  Mach  number  M  =  2.5.  The 
analysis  should  be  repeated  for  the  compressible,  adiabatic,  Falkner- 
Skan  velocity  profiles.  We  have  completed  M  =  0  cases  for  various  values 
of  g,  and  would  like  to  examine  the  influence  of  Mach  number  for  the  same 
values  of  g.  It  was  observed  that  for  g  =  -.0001  and  -.0005  the  conver¬ 
gence  at  the  two  ends  of  the  spectrum  was  very  slow.  These  cases  should 
be  analyzed  somewhat  more  thoroughly. 
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APPENDIX  A 


THE  HOWARD  CIRCLE  THEOREM 


The  Howard  semicircle  theorem  [9]  is  an  extension  of  the  well  known 

fact  that  if  the  amplification  factor  C^  >  0  then  the  propagation  velocity 

C  must  lie  in  the  range  of  U.  Howard  was  able  to  restrict  the  permissible 

values  of  C  and  C^  so  that  the  complex  wave  velocity  C  is  confined  to  a 

semicircle  which  has  the  range  of  U  as  its  diameter.  If  U  and  U  . 

max  mm 

are  the  extrema  of  the  range  of  U,  the  theorem  states 

[Cr  -  l/2(a  +  b)]2  +  C2  <  [l/2(a  +  b)]2  ,  C±  >  0 

where  a  =  U  ,  b  =  U  .  . 

max  min 
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APPENDIX  B 


EIGENVALUES  FROM  STABILITY  ANALYSIS 
FOR  REVERSED  FLOW  BOUNDARY  LAYERS 


EIGENVALUES  FROM  STABILITY  ANALYSIS  FOR 
REVERSED  FLOW  BOUNDARY  LAYERS 

TABLE  B-l 
8  =  -.0001 


a 

C 

r 

C. 

i 

0 

.90538414714 

.025680518247 

.01 

.91223794353 

.071997946349 

.02 

.89422525454 

.11766985541 

.03 

.87074750211 

.12610731927 

.04 

.8412659309 

.14572816584 

.05 

.80533205 

.17071090 

.07 

.76425576 

.21131732 

.10 

.72900063223 

.23958792186 

.15 

.67639882435 

.24398318931 

.18 

.651890026062 

.235453936491 

.20 

.637798788145 

.22696852817 

.22 

.62587353506 

.21686502743 

.25 

.60392754004 

.19700990613 

.27 

.60253175376 

.18702547786 

.29 

.58908572761 

.17144819857 

.30 

.58616022 

.16481124 

.31 

.58992531643 

.16063748234 

.32 

.58748966823 

.15388430488 

.35 

.58168544052 

.13349344451 

.40 

.57632116036 

.099664728940 

.41 

.56994477164 

.090751023415 

.42 

.56963082632391 

.08414875686479 
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EIGENVALUES  FROM  STABILITY  ANALYSIS  FOR 
REVERSED  FLOW  BOUNDARY  LAYERS 


TABLE  B-2 

3  =  -.0005 

Ct 

C 

r 

C. 

i 

0 

.91530377348 

.02166761564 

.01 

,91372171557 

.025240777143 

.05 

.85274628651 

.13020125042 

.10 

.74360436406 

.21320430557 

.15 

.67480093675282 

.24234861164626 

.20 

.63120621598575 

.23023864570417 

.25 

.60282728906963 

.20264037260815 

.30 

.58491658821 

.16939693499 

.35 

.57438769614 

.13455531319 

.40 

.56967589511 

.099924636095 

.45 

.56930353919 

.066316584587 

.46 

.5696486522 

.059749391713 

EIGENVALUES  FROM  STABILITY  ANALYSIS  FOR 
REVERSED  FLOW  BOUNDARY  LAYERS 


a 

TABLE  B-3 

g  =  -.002 

C 

C, 

r 

i 

0 

.9237409069 

.0095488477066 

.01 

.92117830521 

.010983058744 

.05 

.88971797 

.091967207 

.10 

.78686024 

.15214965 

.15 

.70466415 

.20301679 

.20 

.64798726 

.21161601 

.25 

.61069242 

.19445403 

.30 

.58711391 

.16596060 

.35 

.57311136 

.13306074 

.40 

.56593472 

.098968376 

.45 

.56371274 

.065174535 

.50 

.56504181973 

.032280308591 

.55 

.56865534872 

.00023011047136 

.56 

•57775580181 

. 19174876005 (10) _1° 

.57 

.58685371398 

. 61658214056 ( 10) _11 

.58 

.58704040582 

.34648626536U0)"11 
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EIGENVALUES  FROM  STABILITY  ANALYSIS  FOR 
REVERSED  FLOW  BOUNDARY  LAYERS 


a 

TABLE  B-4 

8  =  -.04 

C 

c 

.12 

r 

.9462446953107 

1 

.00077598474441 

.13 

.94480294724055 

.00092398603074 

.14 

.93834475353190 

.0019712426475 

.15 

.92864533 

.0041227836 

.17 

.91398636 

.009481253 

.20 

.88234356 

.031932225 

.23 

.78918067 

.079558300 

.25 

.73325251 

.10055163 

.30 

.63589654 

.12879591 

.35 

.57744885 

.12534423 

.40 

.54449017 

.10193011 

.42 

.53642765 

.089694409 

.44 

.53060095 

.076599668 

.46 

.52663117 

.062983805 

.48 

.52423199 

.049089198 

.50 

.52306536 

.035068680 

.52 

.52292271 

.021022046 

.54 

.52356367 

.0069847632 

.55 

.52416855 

.34183573(10)'6 

.56 

.53445687 

.11840627(10)"8 
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TABLE  Br-4  (con't) 
0  =  -.04 


a 

C 

r 

c. 

1 

57 

.53445745 

.28132369(10)-9 

58 

.53445776 

. 10618506(10)“9 

59 

.53445776 

-. 16007328 (10) ~10 

60 

.534457761505 

- .  581878617 (10)-10 

61 

.53445776292 

-.75694864203(10)“ 
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a 

EIGENVALUES  FROM  STABILITY  ANALYSIS  FOR 

REVERSED  FLOW  BOUNDARY  LAYERS 

TABLE  B-5 

g  =  -.08 

C 

c. 

r 

i 

.20 

.94462051904 

.00060442053448 

.22 

.93168283195 

.0023838799717 

.25 

.91152233835 

.007171373662 

.27 

.89227578 

.014429883 

.30 

.83585489 

.036775313 

.35 

.68705498 

.067058934 

.40 

.57949813 

.076982790 

.45 

.51260920 

.067262598 

.47 

.49694639 

.057479918 
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EIGENVALUES  FROM  STABILITY  ANALYSIS  FOR 
REVERSED  FLOW  BOUNDARY  LAYERS 


TABLE  B-6 

3  =  -.12 


a 

C 

r 

Ci 

.28 

.93771615513 

.00057735106118 

.30 

.92267372951 

.0025385188104 

.32 

.90428630103 

.006241760255 

.35 

.86703422542955 

.016184737774036 

.37 

.82518601942224 

.025407066978775 

.40 

.74110563666309 

.033621469860405 

.42 

.68545236670957 

.034160709745084 

.45 

.60602465101 

.029295418399 

.47 

.54979488217 

.020797858572 

.50 

.47928620449 

.17208372583(10)'' 

.51 

.47900001254 

.164119887 (10)~8 

.52 

.47900001195552 

. 2969542445 (10)"9 
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EIGENVALUES  FROM  STABILITY  ANALYSIS  FOR 
REVERSED  FLOW  BOUNDARY  LAYERS 


TABLE  B-7 

8  =  -.16 

a 

C 

CJ 

r 

i 

35 

.92539848853 

.00064674041646 

37 

.90652445594 

.0030732291515 

40 

.87214263 

.0087996030 

42 

.839073304 

.013608897 

45 

.76943755 

.016917978 

47 

.71714018755 

.014067976705 

50 

.63945116 

.0025996533 

52 

.59255631588 

-.0054550150066 
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EIGENVALUES  FROM  STABILITY  ANALYSIS  FOR 
REVERSED  FLOW  BOUNDARY  LAYERS 


TABLE  B-8 

8  =  -.19884 

a 

C 

r 

Ci 

.37 

.92379057577 

.00067331237442 

.38 

.91876397945 

.00087323704317 

.39 

.91329400663 

.001065562476 

.40 

.90688937 

.0013915674 

.42 

.89032638669 

.00295916137 

.45 

.85880806 

.0055504066 

.47 

.83019514805718 

.0066806913008495 

.50 

.77300900 

.0043585794 

.52 

.7283639558 

.48854083711 (10) 

.53 

.70401641 

.28953714(10)~7 

.54 

.70412941 

-. 19018251 (10) ~7 

r 
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FIGURE  B-l 
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APPENDIX  C 
COMPUTER  PROGRAM 

The  following  FORTRAN  program  was  used  in  the  search  for  eigenvalues. 
The  driver  program  FINDMIN  gives  the  initial  guesses  for  CREAL  and  CIM  and 
then  calls  the  minimization  routine.  The  error  is  returned  as  SSQ.  Sub¬ 
routine  MAINFCN  does  the  integration  using  a  system  routine  RKDF  which 
uses  a  fourth  order  Runge-Kutta  method.  The  arguments  to  RKDF  are: 

X  -  the  independent  variable,  Y  -  the  dependent  variables,  N  -  the  number 
of  variables,  DX  -  the  step  size,  and  IER  -  an  error  return.  RKDF  also 
requires  a  function  F  which  computes  the  derivatives  of  the  dependent 
variables  and  stores  them  in  P. 

The  minimization  routines  are  fairly  general.  The  equivalence  causes 

the  minization  to  be  done  with  respect  to  CREAL  and  CIM.  To  minimize  with 

respect  to  CIM  and  ALPHA  the  equivalence  statement  would  be:  "EQUIVALENCE 

(CIM,  X(l)),  (ALPHA,  X(2))."  Note  that  the  two  variables  which  are 

equivalenced  with  X(l)  and  X(2)  must  be  stored  consecutively  in  memory. 

The  array  Y  represents  the  following  quantities.  Y(l)  =  f,  Y(2)  =  f’, 

Y(3)  =  f",  Y (4)  =  <f>R,  Y (5)  =  (jJj,  Y(6)  =  <j>R',  Y(7)  =  (pj.’ .  Convergence  was 

-5  -9 

generally  achieved  when  minimization  errors  of  10  to  10  occurred  for 
most  cases. 
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PROGRAM  FINOMIN1 INPUT , OUTPUT, TAPE  5  =  1 NPUT, TAPE 2=OUTPUT) 

COMMON  CREAL, CIM, ALPHA, BETA, P30 

COMMON  /92/SSQ  "  -  - - 

CREAL  =  .57194967522 
CIM  =  .0  33942 176531 

ALPHA  =  ,5  - - -  - 

CALL  MINI 

WRITE (2,1) ALPHA, BETA, CREAL, CIM, SSQ 

1  FORMAT ( 50 ( 2H  *>/"  AlFH A  =  “ , El 4 . c , 5 X, "BET A =“, £1 4 . 8 /“  CREAL  =  M,E14.fi,5X  ‘ 

iX,‘'CIM  =  '',E14.8,5X,"EFRQR=‘‘,E14.6/50(2H  *)) 

IOC  CONTINUE 

ENO  . ""  . 

SUBROUTINE  MAINFCN 

COMMON  CREAL, CIM, ALPHA, BETA, P3D 

COMMON  /B2/  SSQ  ~ '  . ~ 

DIMENSION  r (7) , P (7) 

BETA  =  - « 0  30  c 
Y ( 3)  =  - .CO  51 546 
XENO  =  hO. 

X  =  00.0 

Y ( 1 )  =  u.od  '  “ 

Y  <  2)  »  0.0 
Y (4)  =  0.0 

Y  (5)  =0.0 

CCBAR  =  CKEAL*CREAL  ♦  CIM*CIM 
oCC  =  BETA  /  CCflAR 

FACT  R  =  i,L?HA*ALPHA  ♦  30C*CR£AL'  ' 

R4  —  F  ACTR*  F  ACT  R  ♦  BCC* 5CC*CI M*CIM 
R  =  SORT (SORT (R4|) 


gamma  =  0.5*ATAN2<-3CC*CIM,FA'CT  R" ) 

Y (6)  =  R  *  COS(GAMPA) 

Y ( 7)  =  R  ♦  SIN(GAMMA) 

OX  =  0.C5 
N  =  7 

CALL  F(X,Y,P) 

C***  ,*♦**., LOOP  INTEGRATION 

2  CONTINUE 

CALL  RKOF(X,Y,N,DX, IER) 

IF(X.LE.XENO)  GO  TO  2  . . *  . . 

c»**  ,INNER  L00P  INTEGRATION 

SSQ=Y(5)*Y(E)*Y(4) *Y(4) 


RETURN 

ENO 

SUBROUTINE  F(X,Y,P) 

COMMON  CREAL, CIM,  ALFHA ,  BETA',  P30 
DIMENSION  Y(7>  ,  P(7) 

P(l)  =  Y(2) 

P(2)  =  Y  ( 3) 

P  ( 3)  *  BET  A*  (Y  (2 )  -1.IMYC2)  ♦  l.»  -  Y<i>*Y<3» 
P  ( 4)  =  Y  (6) 

P(5)  *  Y(7> 

U  =  Y (2) 

UOB  =  P (3) 

0  *  UOB/<<U-CREAD**2  ♦  CIM**2> 

B  =  0*CIM 

A  =  ALPHA’ALPHA  ♦  D*(U-CREAL) 

P(  6)  =  A  *  Y  (4 )  -  3*Y<5> 

P( 7)  *  B*Y(4)  -  A*Y (2) 

RETURN 

ENO 

SUBROUTINE  NIM1 _ _ 

COMMON  CREAL, CIM, ALPHA, BET  A  ,  P  3  D 


A 
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O  O  o  OOO  O  O  o  O  o  o  OJ  <\l  :  *4000 


COMMON  /32/  SSQ 

DIMENSION  XEST (2,2) ,  X<2> .STEP (2) 

EQUIVALENCS<C^EAl,X (1 )) , (C IM, X ( 2)  » 

STEP  (II  =  l.E-10 
STEP  <  2)  =  1 .  E-<* 

ERS3Q=l.E-3 

P30=  0  .  “ . .  '  ‘  ' 

2  CONTINUE 

DY^STEPli) 

XEST  <1,1>=X<1>  “ 

GRAD  =  1 . E-30 

SEARCH  ALONG  Xl-AXIS"" 

CALL  MIN2(X,DY,GRAO) 

IF(SSQ.LT.ERSSCi)  RETURN . 

XEST ( 1 >  2) =X  < 1) 

ST  EP  < 1 ) =0 Y 
XE  ST  (  2 , 1)  =X  (  2) 

DY=STEP(2) 

GRAO  =  1 .  E«-3C _ _ _ _ _ 

SEARXH  ALONG  X2-AXIS 

CALL  MIN2(X,0Y,GRA0) 

IF(SSQ.LT.ERSSQ) RETURN 
XEST (2,2»=X(2) 

ST  EP ( 2 ) =DY 

GRA0=(XESTtl,2>-XEST(l,l>)/<XEST(2, 2 ) -XEST (2 , 1) ) 

SEARCH  ALONG  LINE 

CALL  M I N2 (X | 3Y  »  GR AO ) 

IF(S3Q.GT.ERSSQ)G0T02 
RETURN 
ENO 

SUBROUTINE  MIN2(X,STEP»GRA0) 

COMMON  /B2/  SSQ 
LOGICAL  OIRPfOIRN 

01 MENSI ON  X <21  ,  Y1 (3) , Y2< 3» , F<3)  . . 

ERSSQ  =  1  •  E-3 

FIND  DIRECTION  ""  ~~  . 

WRITE  <  2 » 20C ) 

0  FORMAT tllX»“FINO  DIRECTION”) 

N=Q 

SGRA0=SIGN(1. jGPAO) 

CALL  MAINFCN 
CONTINUE 
OIRP=. FALSE. 

DIRN=. FALSE. 

X1ST  AR  =  X (1) 

X2ST  AR  =  X(2) _ _ _ _  _ _ _  _ _ 

F<1)*SSQ 
Y 1  ( 1)  =X(1) 

Y2 ( 1) =  X ( 2)  _  _  _  _ _ 

WRITE (2»iOP)  SSQ.Xtl)  ,X(2) 

IF ( SSQ. LT.ERSSQ) RETURN 

100  FORMAT  ( 1QX» ”ERR=_” >  El 7 . 1  l»_"CRf_  "»  E17 • 1JL^“  CJ^“,E17.11> 
CONTINUE 

TRY  POSITIVE  INCREMENT 

ST  EPX  *ST  EP  _ _ 

OX«STEPX/SORT I1.»GRAD**2) 
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XU)  =  X  1ST  AR»DX 

X(2) =  X2STAR«-SGRA0*SQRT (S TEPX** 2-DX* *  2) 
CALL  HAINFCN 
F<2)  =SSQ 
Yl<2) =X(1> 

Y2(2) =X(2) 

IF(F(2>-F(l)>9,ll,ii 

9  CONTINUE 
C 

C  POSITIVE  INCREMENT  WORKED 
C 

DIRPa .TRUE. 

STEPX=2*STEP 

DX  =  ST EPX/SQRT  < 1 • FGRAD* *  2 ) 
X(l)=XlSTARfOX 
X(2)=X2STAR*5GRA0*S0RT(STEPX**2-0X**2) 
CALL  MAINFGN 
F  ( 3) =SSQ 
Y1(3)=XC1)  " 

Y2(3) =  X  (2) 

WRITE(2,10C)  SS0,X(1) ,X(2) 

IF (SSQ.LT.ERSSO) RETURN 
G0T014 

11  CONTINUE 
C  ' 

C  TRY  NEGATIVE  INCREMENT 
C 

STEPX=STEP 

OX=ST EPX/SQRT (l.*GRA0**2) 

X(l) = X 1ST  AR-OX 
X (2) =X2STAR-SGRAO*SORT(STEPX**2-DX**2) 
CALL  HAINFCN 
F  1 2) =SSQ 

Y1(2)=X(1)  . 

Y2(2) =  X ( 2) 

WR1TE(2,10C)  3S3,X(1) ,X<2> 

IF (SSQ.LT.ERSSC) RETURN 
IF(F<2)-F(1)) 10 » 12 » 12 

10  CONTINUE 

OIRN= .TRUE.  . * 

C 

C  NEGATIVE  INCREMENT  WORKEO 
C 

STEPX=2.*STEP 
OX=STEPX/SORT (l.+GRAO**2) 

X(l> =X1STAR-0X 
X(2)=X2STAP-SGR60*SQRT(STlPX**2-DX**2) 
CALL  MAINFCN 
F  <  3) =SSQ 
Y1  (3) =X(1) 

Y2(3) =X ( 2) 

WRITE ( 2 » ICO )  330  >  X  fl)  t  X  < 2)  " 

IF (SSQ.LT.ERSSQ) RETURN 
G0T014  _  _ _ 

12  "  CCNTINUE 
C 

C  NEITHER  WORKEO.  HALVE  STEP 

c  . '  . .  . 

STEP=STEP/2. 

G0T013 

14  CONTINUE 
C 

C  DIRECTION  FOUND 

C 

WRITE  <2*261) 

201  ~  FORMAT  <10X,"3RACKEt  MINIMUM1*) 


33 


o  o  o  o  o 


IF (OlkP) XSIGN: u 
IF  IDIRN) XSIGN=-1 

15  CONTINUE 
IF(F(3)-F<2>)16«17,17 

16  CONTINUE 
N  =  N*1 

ST  E°X  =N»STEp 

DX=ST  EPX/SQRT (1. »GRAD**2) 

X ( 1) =X ( 1) ♦XSIGN’OX 

X ( 2) =X (2) txSIGN«SGRAD*SQRT(STEPX**2-DX**2> 
Y1 £1) =  Y 1  ( 2 ) 
r  1  (2)  =  Y  1  {  3) 

Y1 <  3) =X(1) 

Y  2 ( 1 ) =Y2 (2) 

Y2 (2) =Y2(3) 

Y2 ( 3) =  X  (  2) 

F  ( 1)  =F  (  2) 

F  ( 2)  =  F  (  3) 

CALL  MAINFCN 
F  (  3)  =  SSQ 

WR  IT  E  (  2  >  IC’f  )  SSG, X( 1)  ,  x  m 
IF(SSO.LT.ERSSG) RETURN 
IFIF(3)-F(2>)lb,17,17 

17  CONTINUE 


MINI li'JM  6RACKETTED 
NOW  FIT  QUADRATIC 
WRITE  (2,202) 

202  FCR1AT  (10X,“UEE  QUADRATIC  APPROX  FOR  MINIMUM**) 
IF  <A3S  (GRAD)  .GT.O  .56*10)  GOTO 3 
F 1  =  Y1 ( 1) «Y1 ( 2) 

F2-Y1 (1) -Y1 (3) 

F3=Y1 ( 2 ) - Y 1 (3) 

£Irl  =  F  (D/Fl/FZ 
DIT2--F  (2)/Fi/F3  ' 

TIT3=F (31/F2/F3 

CITl  =  Yl(l)M3IT2*eiT?> 

CIT2  =  Y1  (2)M3IT1*9IT3> 

CIT3=Y1 <3)» <5IT1*PIT2> 

X(l)=(CITl*ClT2*CIT3)/2./(-3ITl+8IT2*BIT3) 

IF  (ASS  (GRAD)  .iT.  1.5  3-10)  GGTOA 
3  CONTINUE 

F1  =  Y2 ( 1) -Y2<2) 

F2  =Y2 ( 1 ) -Y2 ( 3 ) 

F3=Y2(2)-Y2(3) 

BIT1  =  F  (D/F1/F2  _  _ 

8 1 T  2  =  *  F (2) /F1/F3 

BIT3=F(3)/F2/F3 

CIT1=Y2(1)*(8IT2*9IT?) 

CIT2  =  Y212>*  t8ITl*t)IT3> 

CIT3=Y2(3)*(DIT1*BIT2) 

X  <2)  =  <CITl*CI7  2*CIT3)/2./(t3ITl*3IT2*BIT3>  __ 

A  CGNTINUE 

CALL  MAINFCN 

WRITE (2 ,10C ) S3Q, X (1) , X(2>  _  _ _  _ 

IF <S SO. LT.ERSSO) RETURN 
STEP=STEP/2. 

1  CONTINUE  _ _ _ _ 

RETURN 

ENO 
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LIST  OF  SYMBOLS 


C  =  C  +  iC.,  where  C  and  C.  are  real  and  i  =  /  -1 
r  1  r  i 

propagation  velocity 

C.  amplification  factor 

df  U 

f  defined  by  -r-  =  —  ;  dimensionless  velocity  ratio 

n  ue 

m  pressure  gradient  parameter  (eqn  3.10) 

p  pressure 

U  longitudinal  velocity  component 

V  transverse  velocity  component 

a  wave  number 

u  =  a  QZ. 

dn 

3  pressure  gradient  parameter  (eqn  3.10) 

6  boundary  layer  thickness 

* 

o  displacement  thickness 

£  transformed  similarity  variable  (eqn  3.6) 

ri  transformed  similarity  variable  (eqn  3.7) 

<t> (y)  small  perturbation  variable  for  transverse  velocity  (eqn  4.2) 
v  kinematic  viscosity 


Subscripts 

e  external  flow 

x  partial  derivative  with  respect  to  x 

y  partial  derivative  with  respect  to  y 

q  partial  derivative  with  respect  to  r) 

^  small  perturbation  variable  function  of  y 
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